//----------------------------------*-C++-*----------------------------------//
/*!
 * \file   AngularMesh.cc
 * \file   AngularMesh.i.hh
 * \author Jeremy Roberts
 * \date   06/19/2011
 * \brief  AngularMesh member definitions and template instantiations.
 * \note   Copyright (C) 2011 Jeremy Roberts.  Adapted from Denovo.
 */
//---------------------------------------------------------------------------//
// $Rev:: 161                                           $:Rev of last commit
// $Author:: j.alyn.roberts@gmail.com                   $:Author of last commit
// $Date:: 2011-10-24 11:45:46 +0000 (Mon, 24 Oct 2011) $:Date of last commit
//---------------------------------------------------------------------------//

#include <cmath>
#include <iostream>

#include "angle/QuadratureFactory.hh"
#include "angle/Quadrature.hh"
#include "angle/AngularMesh.hh"

#include "utilities/Constants.hh"


namespace slabtran
{

//---------------------------------------------------------------------------//
// CONSTRUCTOR
//---------------------------------------------------------------------------//
/*!
 * \brief Constructor.
 *
 * \param N \f$ S_N\f$ order
 * 
 * \param do_adjoint build octants such that the angles in each octant have
 * the opposite sign in order to do sweeps for adjoint solves (see Octant)
 */
template <class D>
AngularMesh<D>::AngularMesh( const size_type order,
						     int   quad_type,
                             bool  do_adjoint )
    : d_order(order)
    , d_number_octants( std::pow(2,D::dimension) )
    , d_octants( std::pow(2,D::dimension) ) // 2/4/8 for 1D/2D/3D
{
    using util::SP;

    // octant designator
    typedef typename Octant<D>::Octant_Designator Designator;
    
    // make a quadrature creator
    QuadratureFactory<D> creator;

    // set default normalizations
    double norm = 2.0;
    if ( D::dimension > 1 )
    	norm = constants::four_pi;

    typename QuadratureFactory<D>::SP_Quadrature quad =
    		creator.build(quad_type, d_order, norm);
    Check (quad);

    // Build the octants.  Is this the best way? Do I really want a Designator hard coded to dim 3??
    d_octants[PPP] = new Octant<D>(Designator(1, 1, 1), *quad, do_adjoint);
    d_octants[MPP] = new Octant<D>(Designator(-1, 1, 1), *quad, do_adjoint);
    if ( D::dimension > 1 )
    {
    	d_octants[PMP] = new Octant<D>(Designator(1, -1, 1), *quad, do_adjoint);
    	d_octants[MMP] = new Octant<D>(Designator(-1, -1, 1), *quad, do_adjoint);
    }
    if ( D::dimension > 2 )
    {
    	d_octants[PPM] = new Octant<D>(Designator(1, 1, -1), *quad, do_adjoint);
    	d_octants[MPM] = new Octant<D>(Designator(-1, 1, -1), *quad, do_adjoint);
    	d_octants[PMM] = new Octant<D>(Designator(1, -1, -1), *quad, do_adjoint);
    	d_octants[MMM] = new Octant<D>(Designator(-1, -1, -1), *quad, do_adjoint);
    }

}

// instantiations
template class AngularMesh<_1D>;
template class AngularMesh<_2D>;
template class AngularMesh<_3D>;


} // end namespace slabtran

//---------------------------------------------------------------------------//
//                 end of AngularMesh.cc
//---------------------------------------------------------------------------//
